{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# point_defect_static calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "Description updated: 2019-07-26\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The point_defect_static calculation style computes the formation energy of a point defect by comparing the energies of a system before and after a point defect is inserted. The resulting defect system is analyzed using a few different metrics to help characterize if the defect reconfigures to a different structure upon relaxation.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- The computed values of the point defect formation energies and elastic dipole tensors are sensitive to the size of the system.  Larger systems minimize the interaction between the defects, and the affect that the defects have on the system's pressure.  Infinite system formation energies can be estimated by measuring the formation energy for multiple system sizes, and extrapolating to 1/natoms = 0.\n",
    "- Because only a static relaxation is performed, the final configuration might not be the true stable configuration.  Additionally, the stable configuration may not correspond to any of the standard configurations characterized by the point-defect records in the iprPy/library.  Running multiple configurations increases the likelihood of finding the true stable state, but it does not guarantee it.  Alternatively, a dynamic simulation or a genetic algorithm may be more thorough.\n",
    "- The metrics used to identify reconfigurations are not guaranteed to work for all crystals and defects.  Most notably, the metrics assume that the defect's position coincides with a high symmetry site in the lattice.\n",
    "- The current version assumes that the initial defect-free base system is an elemental crystal structure.  The formation energy expression needs to be updated to handle multi-component crystals.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "The method starts with a bulk initial system, and relaxes the atomic positions with a LAMMPS simulation that performs an energy/force minimization.  The cohesive energy, $E_{coh}$, is taken by dividing the system's total energy by the number of atoms in the system.\n",
    "\n",
    "A corresponding defect system is then constructed using the atomman.defect.point() function.  The defect system is relaxed using the same energy/force minimization as was done with the bulk system.  The formation energy of the defect, $E_{f}^{ptd}$, is obtained as\n",
    "\n",
    "$$E_{f}^{ptd} = E_{total}^{ptd} - E_{coh} * N^{ptd},$$\n",
    "\n",
    "where $E_{total}^{ptd}$ is the total potential energy of the relaxed defect system, and $N^{ptd}$ is the number of atoms in the defect system.\n",
    "\n",
    "The elastic dipole tensor, $P_{ij}$, is also estimated for the point defect. $P_{ij}$ is a symmetric second rank tensor that characterizes the elastic nature of the defect.  Here, $P_{ij}$ is estimated using \\[[1](https://doi.org/10.1080/01418618108239410), [2](https://doi.org/10.1080/01418618308244326)\\]\n",
    "\n",
    "$$ P_{ij} = -V \\langle \\sigma_{ij} \\rangle,$$\n",
    "\n",
    "where $V$ is the system cell volume and $\\langle \\sigma_{ij} \\rangle$ is the residual stress on the system due to the defect, which is computed using the measured cell stresses (pressures) of the defect-free system, $\\sigma_{ij}^{0}$, and the defect-containing system, $\\sigma_{ij}^{ptd}$\n",
    "\n",
    "$$\\langle \\sigma_{ij} \\rangle = \\sigma_{ij}^{ptd} - \\sigma_{ij}^{0}.$$\n",
    "\n",
    "The atomman.defect.point() method allows for the generation of four types of point defects:\n",
    "\n",
    "- __vacancy__, where an atom at a specified location is deleted.\n",
    "\n",
    "- __interstitial__, where an extra atom is inserted at a specified location (that does not correspond to an existing atom).\n",
    "\n",
    "- __substitutional__, where the atomic type of an atom at a specified location is changed.\n",
    "\n",
    "- __dumbbell__ interstitial, where an atom at a specified location is replaced by a pair of atoms equidistant from the original atom's position.\n",
    "\n",
    "Point defect complexes (clusters, balanced ionic defects, etc.) can also be constructed piecewise from these basic types.\n",
    "\n",
    "The final defect-containing system is analyzed using a few simple metrics to determine whether or not the point defect configuration has relaxed to a lower energy configuration:\n",
    "\n",
    "- __centrosummation__ adds up the vector positions of atoms relative to the defect's position for all atoms within a specified cutoff. In most simple crystals, the defect positions are associated with high symmetry lattice sites in which the centrosummation about that point in the defect-free system will be zero. If the defect only hydrostatically displaces neighbor atoms, then the centrosummation should also be zero for the defect system. This is computed for all defect types.\n",
    "\n",
    "$$ \\vec{cs} = \\sum_i^N{\\left( \\vec{r}_i - \\vec{r}_{ptd} \\right)} $$\n",
    "\n",
    "- __position_shift__ is the change in position of an interstitial or substitutional atom following relaxation of the system. A non-zero value indicates that the defect atom has moved from its initially ideal position.\n",
    "\n",
    "$$ \\Delta \\vec{r} = \\vec{r}_{ptd} - \\vec{r}_{ptd}^{0}$$\n",
    "\n",
    "- __db_vect_shift__ compares the unit vector associated with the pair of atoms in a dumbbell interstitial before and after relaxation. A non-zero value indicates that the dumbbell has rotated from its ideal configuration.\n",
    "\n",
    "$$ \\Delta \\vec{db} = \\frac{\\vec{r}_{db1} - \\vec{r}_{db2}}{|\\vec{r}_{db1} - \\vec{r}_{db2}|} - \\frac{\\vec{r}_{db1}^0 - \\vec{r}_{db2}^0}{|\\vec{r}_{db1}^0 - \\vec{r}_{db2}^0|}$$\n",
    "\n",
    "If any of the metrics have values not close to (0,0,0), then there was likely an atomic configuration relaxation.\n",
    "\n",
    "The final defect system and the associated perfect base system are retained in the calculation's archive. The calculation's record reports the base system's cohesive energy, the point defect's formation energy, and the values of any of the reconfiguration metrics used.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [DataModelDict](https://github.com/usnistgov/DataModelDict)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2019-07-29 using iprPy version 0.9.0\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import sys\n",
    "import uuid\n",
    "import shutil\n",
    "import datetime\n",
    "from copy import deepcopy\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np \n",
    "\n",
    "# https://github.com/usnistgov/DataModelDict \n",
    "from DataModelDict import DataModelDict as DM\n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'point_defect_static'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ (required) is the LAMMPS command to use.\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential_file__ gives the path to the potential_LAMMPS reference record to use.  Here, this parameter is automatically generated using potential_name and librarydir.\n",
    "\n",
    "- __potential_dir__ gives the path for the folder containing the artifacts associated with the potential (i.e. eam.alloy file).  Here, this parameter is automatically generated using potential_name and librarydir.\n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required).  Here, this parameter is automatically generated from potential_file and potential_dir."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully loaded potential 1999--Mishin-Y--Ni--LAMMPS--ipr1\n"
     ]
    }
   ],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Define potential_file and potential_dir using librarydir and potential_name\n",
    "potential_file = Path(iprPy.libdir, 'potential_LAMMPS', f'{potential_name}.json')\n",
    "potential_dir = Path(iprPy.libdir, 'potential_LAMMPS', potential_name)\n",
    "\n",
    "# Initialize Potential object using potential_file and potential_dir.\n",
    "potential = lmp.Potential(potential_file, potential_dir)\n",
    "print('Successfully loaded potential', potential)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __prototype_name__ gives the name of the crystal_prototype reference record in the iprPy library to load. \n",
    "\n",
    "- __symbols__ is a list of the potential's elemental model symbols to associate with the unique atom types of the loaded system. \n",
    "\n",
    "- __box_parameters__ is a list of the a, b, c lattice constants to assign to the loaded file.\n",
    "\n",
    "- __load_file__ gives the path to the atomic configuration file to load for the ucell system.  Here, this is generated automatically using prototype_name and librarydir.\n",
    "\n",
    "- __load_style__ specifies the format of load_file.  Here, this is automatically set for crystal_prototype records.\n",
    "\n",
    "- __load_options__ specifies any other keyword options for properly loading the load_file.  Here, this is automatically set for crystal_prototype records.\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is generated using the load parameters and symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.520,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.520,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.520]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.760 |   1.760\n",
      "      2 |       1 |   1.760 |   0.000 |   1.760\n",
      "      3 |       1 |   1.760 |   1.760 |   0.000\n"
     ]
    }
   ],
   "source": [
    "prototype_name = 'A1--Cu--fcc'\n",
    "symbols = ['Ni']\n",
    "box_parameters = uc.set_in_units([3.52, 3.52, 3.52], 'angstrom')\n",
    "\n",
    "# Define load_file using librarydir and prototype_name\n",
    "load_file = Path(iprPy.libdir, 'crystal_prototype', f'{prototype_name}.json')\n",
    "\n",
    "# Define load_style and load_options for crystal_prototype records\n",
    "load_style = 'system_model'\n",
    "load_options = {}\n",
    "\n",
    "# Create ucell by loading prototype record\n",
    "ucell = am.load(load_style, load_file, symbols=symbols, **load_options)\n",
    "\n",
    "# Rescale ucell using box_parameters\n",
    "ucell.box_set(a=box_parameters[0], b=box_parameters[1], c=box_parameters[2], scale=True)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Specify the defect parameters\n",
    "\n",
    "- __pointdefect_name__ gives the name of a point-defect reference record in the iprPy library containing point defect input parameters.\n",
    "\n",
    "- __pointdefect_file__ gives the path to a point_defect reference containing point defect input parameters.  Here, this is built automatically using pointdefect_name and librarydir.\n",
    "\n",
    "- __point_kwargs__ (required) is a dictionary or list of dictonaries containing parameters for generating the defect. Here, values are extracted from pointdefect_file. Allowed keywords are:\n",
    "\n",
    "    - __ptd_type__ indicates which defect type to generate: 'v' for vacancy, 'i' for interstitial, 's' for substitutional, or 'db' for dumbbell.\n",
    "    \n",
    "    - __atype__ is the atom type to assign to the defect atom ('i', 's', 'db' ptd_types).\n",
    "    \n",
    "    - __pos__ specifies the position for adding the defect atom (all ptd_types).\n",
    "    \n",
    "    - __ptd_id__ specifies the id of an atom in the initial system where the defect is to be added. Alternative to using pos ('v', 's', 'db' ptd_types).\n",
    "    \n",
    "    - __db_vect__ gives the vector associated with the dumbbell interstitial to generate ('db' ptd_type).\n",
    "    \n",
    "    - __scale__ indicates if pos and db_vect are in absolute (False) or box-relative (True) coordinates. Default is False.\n",
    "    \n",
    "    - __atol__ is the absolute tolerance for position-based searching. Default is 1e-3 angstroms.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "point_kwargs =\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[{'ptd_type': 'v', 'pos': array([0., 0., 0.]), 'scale': False}]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pointdefect_name = 'A1--Cu--fcc--vacancy'\n",
    "#pointdefect_name = 'A1--Cu--fcc--1nn-divacancy'\n",
    "#pointdefect_name = 'A1--Cu--fcc--2nn-divacancy'\n",
    "#pointdefect_name = 'A1--Cu--fcc--100-dumbbell'\n",
    "#pointdefect_name = 'A1--Cu--fcc--110-dumbbell'\n",
    "#pointdefect_name = 'A1--Cu--fcc--111-dumbbell'\n",
    "#pointdefect_name = 'A1--Cu--fcc--octahedral-interstitial'\n",
    "#pointdefect_name = 'A1--Cu--fcc--tetrahedral-interstitial'\n",
    "#pointdefect_name = 'A1--Cu--fcc--crowdion-interstitial'\n",
    "\n",
    "# Define pointdefect_file using librarydir and pointdefect_name\n",
    "pointdefect_file = Path(iprPy.libdir, 'point_defect', f'{pointdefect_name}.json')\n",
    "\n",
    "# Parse pointdefect_file using iprPy.input.interpret()\n",
    "defectinputs = {'ucell':ucell, 'pointdefect_file':pointdefect_file}\n",
    "iprPy.input.subset('pointdefect').interpret(defectinputs)\n",
    "\n",
    "# Extract point_kwargs\n",
    "point_kwargs = defectinputs['point_kwargs']\n",
    "print('point_kwargs =')\n",
    "point_kwargs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Modify system\n",
    "\n",
    "- __sizemults__ list of three integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in creating system.\n",
    "\n",
    "- __system__ is an atomman.System to perform the scan on (required). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# of atoms in system = 4000\n"
     ]
    }
   ],
   "source": [
    "sizemults = [10, 10, 10]\n",
    "\n",
    "# Generate system by supersizing ucell\n",
    "system = ucell.supersize(*sizemults)\n",
    "print('# of atoms in system =', system.natoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.6. Specify calculation-specific run parameters\n",
    "\n",
    "- __energytolerance__ is the energy tolerance to use during the minimizations. This is unitless.\n",
    "\n",
    "- __forcetolerance__ is the force tolerance to use during the minimizations. This is in energy/length units.\n",
    "\n",
    "- __maxiterations__ is the maximum number of minimization iterations to use.\n",
    "\n",
    "- __maxevaluations__ is the maximum number of minimization evaluations to use.\n",
    "\n",
    "- __maxatommotion__ is the largest distance that an atom is allowed to move during a minimization iteration. This is in length units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "energytolerance = 1e-8\n",
    "forcetolerance = uc.set_in_units(0.0, 'eV/angstrom')\n",
    "maxiterations = 10000\n",
    "maxevaluations = 100000\n",
    "maxatommotion = uc.set_in_units(0.01, 'angstrom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. min.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('min.template', 'w') as f:\n",
    "    f.write(\"\"\"# LAMMPS input script that performs a simple energy minimization\n",
    "\n",
    "box tilt large\n",
    "\n",
    "<atomman_system_info>\n",
    "\n",
    "<atomman_pair_info>\n",
    "\n",
    "thermo_style custom step lx ly lz pxx pyy pzz pxy pxz pyz pe\n",
    "thermo_modify format float %.13e\n",
    "\n",
    "compute peatom all pe/atom \n",
    "\n",
    "dump dumpit all custom <maxeval> atom.* id type x y z c_peatom\n",
    "dump_modify dumpit format <dump_modify_format>\n",
    "\n",
    "min_modify dmax <dmax>\n",
    "\n",
    "minimize <etol> <ftol> <maxiter> <maxeval>\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2. pointdefect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pointdefect(lammps_command, system, potential, point_kwargs,\n",
    "                mpi_command=None, etol=0.0, ftol=0.0, maxiter=10000,\n",
    "                maxeval=100000, dmax=uc.set_in_units(0.01, 'angstrom')):\n",
    "    \"\"\"\n",
    "    Adds one or more point defects to a system and evaluates the defect \n",
    "    formation energy.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    point_kwargs : dict or list of dict\n",
    "        One or more dictionaries containing the keyword arguments for\n",
    "        the atomman.defect.point() function to generate specific point\n",
    "        defect configuration(s).\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    sim_directory : str, optional\n",
    "        The path to the directory to perform the simuation in.  If not\n",
    "        given, will use the current working directory.\n",
    "    etol : float, optional\n",
    "        The energy tolerance for the structure minimization. This value is\n",
    "        unitless. (Default is 0.0).\n",
    "    ftol : float, optional\n",
    "        The force tolerance for the structure minimization. This value is in\n",
    "        units of force. (Default is 0.0).\n",
    "    maxiter : int, optional\n",
    "        The maximum number of minimization iterations to use (default is \n",
    "        10000).\n",
    "    maxeval : int, optional\n",
    "        The maximum number of minimization evaluations to use (default is \n",
    "        100000).\n",
    "    dmax : float, optional\n",
    "        The maximum distance in length units that any atom is allowed to relax\n",
    "        in any direction during a single minimization iteration (default is\n",
    "        0.01 Angstroms).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'E_coh'** (*float*) - The cohesive energy of the bulk system.\n",
    "        - **'E_ptd_f'** (*float*) - The point.defect formation energy.\n",
    "        - **'E_total_base'** (*float*) - The total potential energy of the\n",
    "          relaxed bulk system.\n",
    "        - **'E_total_ptd'** (*float*) - The total potential energy of the\n",
    "          relaxed defect system.\n",
    "        - **'system_base'** (*atomman.System*) - The relaxed bulk system.\n",
    "        - **'system_ptd'** (*atomman.System*) - The relaxed defect system.\n",
    "        - **'dumpfile_base'** (*str*) - The filename of the LAMMPS dump file\n",
    "          for the relaxed bulk system.\n",
    "        - **'dumpfile_ptd'** (*str*) - The filename of the LAMMPS dump file\n",
    "          for the relaxed defect system.\n",
    "    \"\"\"\n",
    "    try:\n",
    "        # Get script's location if __file__ exists\n",
    "        script_dir = Path(__file__).parent\n",
    "    except:\n",
    "        # Use cwd otherwise\n",
    "        script_dir = Path.cwd()\n",
    "\n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "    \n",
    "    #Get lammps version date\n",
    "    lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "    \n",
    "    # Define lammps variables\n",
    "    lammps_variables = {}\n",
    "    system_info = system.dump('atom_data', f='perfect.dat',\n",
    "                              units=potential.units,\n",
    "                              atom_style=potential.atom_style)\n",
    "    lammps_variables['atomman_system_info'] = system_info\n",
    "    lammps_variables['atomman_pair_info'] = potential.pair_info(system.symbols)\n",
    "    lammps_variables['etol'] = etol\n",
    "    lammps_variables['ftol'] = uc.get_in_units(ftol, lammps_units['force'])\n",
    "    lammps_variables['maxiter'] = maxiter\n",
    "    lammps_variables['maxeval'] = maxeval\n",
    "    lammps_variables['dmax'] = dmax\n",
    "    \n",
    "    # Set dump_modify_format based on lammps_date\n",
    "    if lammps_date < datetime.date(2016, 8, 3):\n",
    "        lammps_variables['dump_modify_format'] = '\"%d %d %.13e %.13e %.13e %.13e %.13e %.13e %.13e\"'\n",
    "    else:\n",
    "        lammps_variables['dump_modify_format'] = 'float %.13e'\n",
    "    \n",
    "    # Write lammps input script\n",
    "    template_file = Path(script_dir, 'min.template')\n",
    "    lammps_script = 'min.in'\n",
    "    with open(template_file) as f:\n",
    "        template = f.read()\n",
    "    with open(lammps_script, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables, '<', '>'))\n",
    "\n",
    "    # Run lammps to relax perfect.dat\n",
    "    output = lmp.run(lammps_command, lammps_script, mpi_command)\n",
    "    \n",
    "    # Extract LAMMPS thermo data.\n",
    "    thermo = output.simulations[0]['thermo']\n",
    "    E_total_base = uc.set_in_units(thermo.PotEng.values[-1],\n",
    "                                   lammps_units['energy'])\n",
    "    E_coh = E_total_base / system.natoms\n",
    "    \n",
    "    pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])\n",
    "    pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])\n",
    "    pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])\n",
    "    pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])\n",
    "    pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])\n",
    "    pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])\n",
    "    pressure_base = np.array([[pxx, pxy, pxz], [pxy, pyy, pyz], [pxz, pyz, pzz]])\n",
    "    \n",
    "    # Rename log file\n",
    "    shutil.move('log.lammps', 'min-perfect-log.lammps')\n",
    "    \n",
    "    # Load relaxed system from dump file and copy old box vectors because \n",
    "    # dump files crop the values.\n",
    "    last_dump_file = 'atom.' + str(thermo.Step.values[-1])\n",
    "    system_base = am.load('atom_dump', last_dump_file, symbols=system.symbols)\n",
    "    system_base.box_set(vects=system.box.vects)\n",
    "    system_base.dump('atom_dump', f='perfect.dump')\n",
    "    \n",
    "    # Add defect(s)\n",
    "    system_ptd = deepcopy(system_base)\n",
    "    if not isinstance(point_kwargs, (list, tuple)):\n",
    "        point_kwargs = [point_kwargs]\n",
    "    for pkwargs in point_kwargs:\n",
    "        system_ptd = am.defect.point(system_ptd, **pkwargs)\n",
    "    \n",
    "    # Update lammps variables\n",
    "    system_info = system_ptd.dump('atom_data', f='defect.dat',\n",
    "                                  units = potential.units, \n",
    "                                  atom_style = potential.atom_style)\n",
    "    lammps_variables['atomman_system_info'] = system_info\n",
    "    \n",
    "    # Write lammps input script\n",
    "    with open(lammps_script, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables,\n",
    "                                         '<', '>'))\n",
    "    \n",
    "    # Run lammps\n",
    "    output = lmp.run(lammps_command, lammps_script, mpi_command)\n",
    "    \n",
    "    # Extract lammps thermo data\n",
    "    thermo = output.simulations[0]['thermo']\n",
    "    E_total_ptd = uc.set_in_units(thermo.PotEng.values[-1],\n",
    "                                  lammps_units['energy'])\n",
    "    pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])\n",
    "    pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])\n",
    "    pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])\n",
    "    pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])\n",
    "    pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])\n",
    "    pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])\n",
    "    pressure_ptd = np.array([[pxx, pxy, pxz], [pxy, pyy, pyz], [pxz, pyz, pzz]])\n",
    "    \n",
    "    # Rename log file\n",
    "    shutil.move('log.lammps', 'min-defect-log.lammps')\n",
    "    \n",
    "    # Load relaxed system from dump file and copy old vects as \n",
    "    # the dump files crop the values\n",
    "    last_dump_file = 'atom.'+str(thermo.Step.values[-1])\n",
    "    system_ptd = am.load('atom_dump', last_dump_file, symbols=system_ptd.symbols)\n",
    "    system_ptd.box_set(vects=system.box.vects)\n",
    "    system_ptd.dump('atom_dump', f='defect.dump')\n",
    "    \n",
    "    # Compute defect formation energy\n",
    "    E_ptd_f = E_total_ptd - E_coh * system_ptd.natoms\n",
    "    \n",
    "    # Compute strain tensor\n",
    "    pij = -(pressure_base - pressure_ptd) * system_base.box.volume\n",
    "    \n",
    "    # Cleanup files\n",
    "    for fname in script_dir.glob('atom.*'):\n",
    "        fname.unlink()\n",
    "    for dumpjsonfile in script_dir.glob('*.dump.json'):\n",
    "        dumpjsonfile.unlink()\n",
    "    \n",
    "    # Return results\n",
    "    results_dict = {}\n",
    "    results_dict['E_coh'] = E_coh\n",
    "    results_dict['E_ptd_f'] = E_ptd_f\n",
    "    results_dict['E_total_base'] = E_total_base\n",
    "    results_dict['E_total_ptd'] = E_total_ptd\n",
    "    results_dict['pij_tensor'] = pij\n",
    "    results_dict['system_base'] = system_base\n",
    "    results_dict['system_ptd'] = system_ptd\n",
    "    results_dict['dumpfile_base'] = 'perfect.dump'\n",
    "    results_dict['dumpfile_ptd'] = 'defect.dump'\n",
    "    \n",
    "    return results_dict\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3. check_ptd_config()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_ptd_config(system, point_kwargs, cutoff,\n",
    "                     tol=uc.set_in_units(1e-5, 'angstrom')):\n",
    "    \"\"\"\n",
    "    Evaluates a relaxed system containing a point defect to determine if the\n",
    "    defect structure has transformed to a different configuration.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    system : atomman.System\n",
    "        The relaxed defect system.\n",
    "    point_kwargs : dict or list of dict\n",
    "        One or more dictionaries containing the keyword arguments for\n",
    "        the atomman.defect.point() function to generate specific point\n",
    "        defect configuration(s).\n",
    "    cutoff : float\n",
    "        Cutoff distance to use in identifying neighbor atoms.\n",
    "    tol : float, optional\n",
    "        Absolute tolerance to use for identifying if a defect has\n",
    "        reconfigured (default is 1e-5 Angstoms).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'has_reconfigured'** (*bool*) - Flag indicating if the structure\n",
    "          has been identified as relaxing to a different defect configuration.\n",
    "        - **'centrosummation'** (*numpy.array of float*) - The centrosummation\n",
    "          parameter used for evaluating if the configuration has relaxed.\n",
    "        - **'position_shift'** (*numpy.array of float*) - The position_shift\n",
    "          parameter used for evaluating if the configuration has relaxed.\n",
    "          Only given for interstitial and substitutional-style defects.\n",
    "        - **'db_vect_shift'** (*numpy.array of float*) - The db_vect_shift\n",
    "          parameter used for evaluating if the configuration has relaxed.\n",
    "          Only given for dumbbell-style defects.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Check if point_kwargs is a list\n",
    "    if not isinstance(point_kwargs, (list, tuple)):\n",
    "        pos = point_kwargs['pos']\n",
    "    \n",
    "    # If it is a list of 1, use that set\n",
    "    elif len(point_kwargs) == 1:\n",
    "        point_kwargs = point_kwargs[0]\n",
    "        pos = point_kwargs['pos']\n",
    "        \n",
    "    # If it is a list of two (divacancy), use the first and average position\n",
    "    elif len(point_kwargs) == 2:\n",
    "        pos = (np.array(point_kwargs[0]['pos'])\n",
    "               + np.array(point_kwargs[1]['pos'])) / 2\n",
    "        point_kwargs = point_kwargs[0]\n",
    "    \n",
    "    # More than two not supported by this function\n",
    "    else:\n",
    "        raise ValueError('Invalid point defect parameters')\n",
    "\n",
    "    # Initially set has_reconfigured to False\n",
    "    has_reconfigured = False\n",
    "    \n",
    "    # Calculate distance of all atoms from defect position\n",
    "    pos_vects = system.dvect(system.atoms.pos, pos) \n",
    "    pos_mags = np.linalg.norm(pos_vects, axis=1)\n",
    "    \n",
    "    # Calculate centrosummation by summing up the positions of the close atoms\n",
    "    centrosummation = np.sum(pos_vects[pos_mags < cutoff], axis=0)\n",
    "    \n",
    "    if not np.allclose(centrosummation, np.zeros(3), atol=tol):\n",
    "        has_reconfigured = True\n",
    "        \n",
    "    # Calculate shift of defect atom's position if interstitial or substitutional\n",
    "    if point_kwargs['ptd_type'] == 'i' or point_kwargs['ptd_type'] == 's':\n",
    "        position_shift = system.dvect(system.natoms-1, pos)\n",
    "       \n",
    "        if not np.allclose(position_shift, np.zeros(3), atol=tol):\n",
    "            has_reconfigured = True\n",
    "        \n",
    "        return {'has_reconfigured': has_reconfigured,\n",
    "                'centrosummation': centrosummation,\n",
    "                'position_shift': position_shift}\n",
    "        \n",
    "    # Investigate if dumbbell vector has shifted direction \n",
    "    elif point_kwargs['ptd_type'] == 'db':\n",
    "        db_vect = point_kwargs['db_vect'] / np.linalg.norm(point_kwargs['db_vect'])\n",
    "        new_db_vect = system.dvect(-2, -1)\n",
    "        new_db_vect = new_db_vect / np.linalg.norm(new_db_vect)\n",
    "        db_vect_shift = db_vect - new_db_vect\n",
    "        \n",
    "        if not np.allclose(db_vect_shift, np.zeros(3), atol=tol):\n",
    "            has_reconfigured = True\n",
    "        \n",
    "        return {'has_reconfigured': has_reconfigured,\n",
    "                'centrosummation': centrosummation,\n",
    "                'db_vect_shift': db_vect_shift}\n",
    "    \n",
    "    else:\n",
    "        return {'has_reconfigured': has_reconfigured,\n",
    "                'centrosummation': centrosummation}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1. Generate point defect system and evaluate the energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = pointdefect(lammps_command, system, potential, point_kwargs,\n",
    "                           mpi_command = mpi_command,\n",
    "                           etol = energytolerance,\n",
    "                           ftol = forcetolerance,\n",
    "                           maxiter = maxiterations,\n",
    "                           maxeval = maxevaluations,\n",
    "                           dmax = maxatommotion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['E_coh', 'E_ptd_f', 'E_total_base', 'E_total_ptd', 'pij_tensor', 'system_base', 'system_ptd', 'dumpfile_base', 'dumpfile_ptd'])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.2. Characterize if the defect has reconfigured"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoff = 1.05*ucell.box.a\n",
    "results_dict.update(check_ptd_config(results_dict['system_ptd'], \n",
    "                                     point_kwargs, \n",
    "                                     cutoff))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['E_coh', 'E_ptd_f', 'E_total_base', 'E_total_ptd', 'pij_tensor', 'system_base', 'system_ptd', 'dumpfile_base', 'dumpfile_ptd', 'has_reconfigured', 'centrosummation'])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1. Define units for outputting values\n",
    "\n",
    "- __length_unit__ is the unit of length to display relaxed lattice constants in.\n",
    "- __energy_unit__ is the unit of energy to display cohesive energies in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "length_unit = 'angstrom'\n",
    "energy_unit = 'eV'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2. Print $E_{coh}$ and $E_{ptd}^f$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ecoh =   -4.449999998346 eV\n",
      "Eptd_f = 1.5999212636525046 eV\n"
     ]
    }
   ],
   "source": [
    "print('Ecoh =  ', uc.get_in_units(results_dict['E_coh'], energy_unit), energy_unit)\n",
    "print('Eptd_f =', uc.get_in_units(results_dict['E_ptd_f'], energy_unit), energy_unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3. Check configuration parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Has the system (likely) reconfigured? False\n",
      "centrosummation = [-1.89367102e-12 -1.89381844e-12 -1.89359639e-12] angstrom\n"
     ]
    }
   ],
   "source": [
    "print('Has the system (likely) reconfigured?', results_dict['has_reconfigured'])\n",
    "if 'centrosummation' in results_dict:\n",
    "    print('centrosummation =', uc.get_in_units(results_dict['centrosummation'], length_unit), length_unit)\n",
    "if 'position_shift' in results_dict:\n",
    "    print('position_shift = ', uc.get_in_units(results_dict['position_shift'], length_unit), length_unit)\n",
    "if 'db_vect_shift' in results_dict:\n",
    "    print('db_vect_shift =  ', uc.get_in_units(results_dict['db_vect_shift'], length_unit), length_unit)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
